#ifndef SD_SPARSE_HPP
#define SD_SPARSE_HPP

// -----------------------------------------------------------------------------------------
// Header files
// -----------------------------------------------------------------------------------------
#include "config.hpp"
#include <algorithm>

namespace MT {

//! generate a sparser weight distribution by removing a percentage of small indices
//! Input:
//!   x: a sparse vector representing the current weight distribution
//!   alpha: [0, 1]: percentage to keep
//! Output:
//!   another sparse vector
//! TODO: use argsort() instead and then UnitTest
template<typename _Scalar>
Eigen::SparseVector<_Scalar> sparser(Eigen::SparseVector<_Scalar> const &x, double percentage=0.95)
{
    // get the weights
    int i, j, N = x.m_data.size();
    std::vector<_Scalar> W(N);
    for(i = 0; i < N; ++i) W[i] = x.m_data.value(i);

    // sort them
    std::sort(W.begin(), W.end());

    // find threshold
    _Scalar thresh, w, sw = x.sum() * (1-percentage); // total weight
    for(i = 0, w = 0;  i < N && w <= sw; w += W[i++]);
    thresh = W[--i];

    // count the number of elements greater than thresh
    for(i = 0, j = 0; i < N; ++i) if(x.m_data.value(i) >= thresh) ++j;

    // create result
    Eigen::SparseVector<_Scalar> y(x.size());
    y.m_data.resize(j);
    for(i = 0, j = 0; i < N; ++i) if(x.m_data.value(i) >= thresh)
    {
        y.m_data.index(j) = x.m_data.index(i);
        y.m_data.value(j++) = x.m_data.value(i);
    }

    return y;
}

} // namespace MT

#endif // SD_WEIGHTED_INDEX_HPP
